Apparatus and method for providing near-optimal representations over redundant dictionaries

ABSTRACT

Certain exemplary embodiments provide a method comprising a plurality of activities, comprising: automatically: receiving a signal; and via a plurality of iterations, determining, for the received signal, a provable near-optimal B-term representation formed from a linear combination of a plurality of elements selected from a predetermined redundant μ-coherent dictionary, said determined near-optimal B-term representation corresponding to an error that is a minimum of a set of errors corresponding to all possible B-term representations associated with the predetermined redundant μ-coherent dictionary.

CROSS-REFERENCES TO RELATED APPLICATIONS

[0001] This application claims priority to, and incorporates by reference herein in its entirety, pending U.S. Provisional Patent Application Ser. No. 60/472,614 (Attorney Docket No. 2002-0139), filed 22 May 2003.

SUMMARY

[0002] Certain exemplary embodiments provide a method comprising a plurality of activities, comprising: automatically: receiving a signal; and via a plurality of iterations, determining, for the received signal, a provable near-optimal B-term representation formed from a linear combination of a plurality of elements selected from a predetermined redundant μ-coherent dictionary, said determined near-optimal B-term representation corresponding to an error that is a minimum of a set of errors corresponding to all possible B-term representations associated with the predetermined redundant μ-coherent dictionary.

BRIEF DESCRIPTION OF THE DRAWINGS

[0003] A wide variety of potential embodiments will be more readily understood through the following detailed description, with reference to the accompanying drawings in which:

[0004]FIG. 1 is a block diagram of an exemplary embodiment of a telecommunications system 1000;

[0005]FIG. 2 is a flow diagram of an exemplary embodiment of a method 2000; and

[0006]FIG. 3 is a block diagram of an exemplary embodiment of an information device 3000.

Definitions

[0007] When the following terms are used herein, the accompanying definitions apply:

[0008] automatically—acting or operating in a manner essentially independent of external influence or control. For example, an automatic light switch can turn on upon “seeing” a person in its view, without the person manually operating the light switch.

[0009] comprising—including but not limited to.

[0010] firmware—machine-readable instructions that are stored in a read-only memory (ROM). ROM's can comprise PROMs and EPROMs.

[0011] haptic—both the human sense of kinesthetic movement and the human sense of touch. Among the many potential haptic experiences are numerous sensations, body-positional differences in sensations, and time-based changes in sensations that are perceived at least partially in non-visual, non-audible, and non-olfactory manners, including the experiences of tactile touch (being touched), active touch, grasping, pressure, friction, traction, slip, stretch, force, torque, impact, puncture, vibration, motion, acceleration, jerk, pulse, orientation, limb position, gravity, texture, gap, recess, viscosity, pain, itch, moisture, temperature, thermal conductivity, and thermal capacity.

[0012] I/O device—any sensory-oriented input and/or output device, such as an audio, visual, haptic, olfactory, and/or taste-oriented device, including, for example, a monitor, display, projector, overhead display, keyboard, keypad, mouse, trackball, joystick, gamepad, wheel, touchpad, touch panel, pointing device, microphone, speaker, video camera, camera, scanner, printer, haptic device, vibrator, tactile simulator, and/or tactile pad, potentially including a port to which an I/O device can be attached or connected.

[0013] information device—any device capable of processing information, such as any general purpose and/or special purpose computer, such as a personal computer, workstation, server, minicomputer, mainframe, supercomputer, computer terminal, laptop, wearable computer, and/or Personal Digital Assistant (PDA), mobile terminal, Bluetooth device, communicator, “smart” phone (such as a Handspring Treo-like device), messaging service (e.g., Blackberry) receiver, pager, facsimile, cellular telephone, a traditional telephone, telephonic device, a programmed microprocessor or microcontroller and/or peripheral integrated circuit elements, an ASIC or other integrated circuit, a hardware electronic logic circuit such as a discrete element circuit, and/or a programmable logic device such as a PLD, PLA, FPGA, or PAL, or the like, etc. In general any device on which resides a finite state machine capable of implementing at least a portion of a method, structure, and/or or graphical user interface described herein may be used as an information device. An information device can include well-known components such as one or more network interfaces, one or more processors, one or more memories containing instructions, and/or one or more input/output (I/O) devices, one or more user interfaces, etc.

[0014] Internet—an interconnected global collection of networks that connect information devices.

[0015] iteration—a repetition.

[0016] machine-readable medium—a physical structure from which a machine can obtain data and/or information. Examples include memory devices, punch cards, etc.

[0017] memory device—any device capable of storing analog or digital information, for example, a non-volatile memory, volatile memory, Random Access Memory, RAM, Read Only Memory, ROM, flash memory, magnetic media, a hard disk, a floppy disk, a magnetic tape, an optical media, an optical disk, a compact disk, a CD, a digital versatile disk, a DVD, and/or a raid array, etc. The memory device can be coupled to a processor and can store instructions adapted to be executed by the processor according to an embodiment disclosed herein.

[0018] near-optimal—an approximation that minimizes an error.

[0019] network—a communicatively coupled plurality of communication devices.

[0020] network interface—any device, system, or subsystem capable of coupling an information device to a network. For example, a network interface can be a telephone, cellular phone, cellular modem, telephone data modem, fax modem, wireless transceiver, ethernet card, cable modem, digital subscriber line interface, bridge, hub, router, or other similar device.

[0021] processor—a device for processing machine-readable instruction. A processor can be a central processing unit, a local processor, a remote processor, parallel processors, and/or distributed processors, etc. The processor can be a general-purpose microprocessor, such the Pentium III series of microprocessors manufactured by the Intel Corporation of Santa Clara, Calif. In another embodiment, the processor can be an Application Specific Integrated Circuit (ASIC) or a Field Programmable Gate Array (FPGA) that has been designed to implement in its hardware and/or firmware at least a part of an embodiment disclosed herein.

[0022] system—a collection of devices and/or instructions, the collection designed to perform one or more specific functions.

[0023] user interface—any device for rendering information to a user and/or requesting information from the user. A user interface includes at least one of textual, graphical, audio, video, animation, and/or haptic elements. A textual element can be provided, for example, by a printer, monitor, display, projector, etc. A graphical element can be provided, for example, via a monitor, display, projector, and/or visual indication device, such as a light, flag, beacon, etc. An audio element can be provided, for example, via a speaker, microphone, and/or other sound generating and/or receiving device. A video element or animation element can be provided, for example, via a monitor, display, projector, and/or other visual device. A haptic element can be provided, for example, via a very low frequency speaker, vibrator, tactile stimulator, tactile pad, simulator, keyboard, keypad, mouse, trackball, joystick, gamepad, wheel, touchpad, touch panel, pointing device, and/or other haptic device, etc. A user interface can include one or more textual elements such as, for example, one or more letters, number, symbols, etc. A user interface can include one or more graphical elements such as, for example, an image, photograph, drawing, icon, window, title bar, panel, sheet, tab, drawer, matrix, table, form, calendar, outline view, frame, dialog box, static text, text box, list, pick list, pop-up list, pull-down list, menu, tool bar, dock, check box, radio button, hyperlink, browser, button, control, palette, preview panel, color wheel, dial, slider, scroll bar, cursor, status bar, stepper, and/or progress indicator, etc. A textual and/or graphical element can be used for selecting, programming, adjusting, changing, specifying, etc. an appearance, background color, background style, border style, border thickness, foreground color, font, font style, font size, alignment, line spacing, indent, maximum data length, validation, query, cursor type, pointer type, autosizing, position, and/or dimension, etc. A user interface can include one or more audio elements such as, for example, a volume control, pitch control, speed control, voice selector, and/or one or more elements for controlling audio play, speed, pause, fast forward, reverse, etc. A user interface can include one or more video elements such as, for example, elements controlling video play, speed, pause, fast forward, reverse, zoom-in, zoom-out, rotate, and/or tilt, etc. A user interface can include one or more animation elements such as, for example, elements controlling animation play, pause, fast forward, reverse, zoom-in, zoom-out, rotate, tilt, color, intensity, speed, frequency, appearance, etc. A user interface can include one or more haptic elements such as, for example, elements utilizing tactile stimulus, force, pressure, vibration, motion, displacement, temperature, etc.

[0024] wireless—any means to transmit a signal that does not require the use of a wire or guide connecting a transmitter and a receiver, such as radio waves, electromagnetic signals at any frequency, lasers, microwaves, etc., but excluding purely visual signaling, such as semaphore, smoke signals, sign language, etc.

[0025] wireline—any means to transmit a signal comprising the use of a wire or waveguide (e.g., optical fiber) connecting a transmitter and receiver. Wireline communications can comprise, for example, telephone communications over a POTS network.

DETAILED DESCRIPTION

[0026] 1. Introduction

[0027] One of the central problems of modern mathematical approximation theory is to approximate functions, or signals, concisely, with elements from a large candidate set called a dictionary. Formally, we are given a signal A ε R^(N) and a dictionary D={φ_(i)}_(iε1) is, of unit vectors that span R^(N). A representation R of B terms for input A ε R^(N) is a linear combination of dictionary elements, R=Σ_(iεΛ)α_(i)φ_(i) for φ_(t) ε D and some Λ, |Λ|≦B. Typically, B<<N, so that R is a concise approximation to signal A. The error of the representation indicates by how well it approximates A, and is given by ∥A−R∥₂={square root}{square root over (Σ_(t)|A[t]−R[t]|²)}. The problem is to find the best B-term representation, i.e., find a R that minimizes ∥A−R∥₂ dictionary may be redundant in the sense that there is more than one possible exact representation for A, i.e., |D|>N=dim(R^(N)). Redundant dictionaries are used because, both theoretically and in practice, for important classes of signals, as the size of a dictionary increases, the error and the conciseness of the approximations improve.

[0028] We present the first known efficient algorithm for finding a provably approximate representation for an input signal over redundant dictionaries. We identify and focus on redundant dictionaries with small coherence (i.e., vectors are nearly orthogonal). We present an “algorithm that preprocesses any such dictionary in time and space polynomial in |D|, and obtains an 1+ε approximate representation of the given signal in time nearly linear in signal size N and polylogarithmic in |D |; by contrast, most algorithms in the literature require Ω(|D|) time, and, yet, provide no provable bounds. The technical crux of our result is our proof that two commonly used local search techniques, when combined appropriately, gives a provably near-optimal signal representation over redundant dictionaries with small coherence. Our result immediately applies to several specific redundant dictionaries considered by the domain experts thus far.” In addition, we present new redundant dictionaries which have small coherence (and therefore are amenable to our algorithms) and yet have significantly large sizes, thereby adding to the redundant dictionary construction literature.

[0029] Work with redundant dictionaries forms the emerging field of highly nonlinear approximation theory. We have presented algorithmic results for some of the most basic problems in this area, but other mathematical and algorithmic questions remain to be explored.

[0030] 1.2 Mathematical Approximation Theory

[0031] The main problem of mathematical approximation theory is to approximate functions compactly, i.e., in small space, or using a “small number of terms.” Formally, we are given a signal A ε R^(N) and a dictionary D={φ_(i)}_(iε1) of unit vectors that span R^(N). A representation R of B terms for input A ε R^(N) is a linear combination of dictionary elements, R=Σ_(iεΛ)α_(i)φ_(i) for φ_(i) ε D and some Λ, |Λ|≦B. Typically, B<<N. The error of the representation indicates how well it approximates A. Following the most common practice, we use l² norm: the error of approximation is ∥A−R∥₂={square root}{square root over (Σ_(t)|A[t]−R[t]|²)} henceforth written ∥A−R∥ or (equivalently) its square.

[0032] Mathematical approximation theory has applications in numerical computations, e.g., adaptive partial differential equation solvers, audio signal compression, image compression and statistical estimation with its applications to classification. It is a rich area with significant mathematical achievements and successful applications, both classical and novel; see References 8 and 9 for good surveys, and The Journal of Approximation Theory (Reference 16) and Constructive Approximation (Reference 2) for current mathematical developments.

[0033] There are two approaches to mathematical approximation theory:

[0034] 1. Linear Approximation Theory. We approximate functions using a fixed linear subspace of the dictionary independent of the signal. For example, if the dictionary is the Fourier basis, a B-term approximation of A is given by the first (lowest) B frequencies of its Fourier expansion.

[0035] 2. Nonlinear Approximation Theory. We seek the best or optimal B-term approximation, i.e., R with |Λ|=B such that ∥A−R∥ is minimized. In this setting, the terms used depend on the signal and do not come from a fixed subspace. There are two further variations. In standard nonlinear approximation theory, the dictionary is an orthonormal basis of size N and each function has a unique exact representation R. It is easy to prove by Parseval's equality that the best B term representation comprises the B largest |

A, φ

| over φ ε D. In highly nonlinear approximation theory, the dictionary is redundant and is larger than N. Signals have more than one exact representation over D and the best collection of B terms need not be the B largest.

[0036] Generally, one wants to relate the quality of a representation to other parameters, such as B, the number of terms in the representation (a central mathematical question) or computational cost (in terms of N, B and |D|, a central algorithmic question). There is considerable mathematical analysis of linear and standard nonlinear approximation theory and the algorithmic issues in both are fairly straightforward. Generalizing from linear to standard non-linear to highly-non-linear theory yields much better representations, but pose much harder mathematical and algorithmic challenges. “Redundancy on the one hand offers much promise for greater efficiency in terms of approximation rate, but on the other hand, gives rise to highly nontrivial theoretical and practical problems” (quoted from Reference 23). In fact, the mathematics of the theory of highly nonlinear approximation is only now emerging, and algorithmic complexity of these problems is wide open. The subject of this paper is the central algorithmic issue in highly non-linear approximation theory: for a given function A, parameter B and a prespecified dictionary D, finding the optimal representation of B terms.

[0037] The general highly-non-linear approximation problem of our interest is NP-hard (see Reference 7). In fact, the proof there implicitly shows that the problem is NP-hard even to determine whether the optimal error is zero for the given B; hence, unless P=NP, no polynomial time algorithm exists that approximates the best B-term representation over an arbitrary dictionary even if we wish only to approximate the optimal error by a factor, however large.

[0038] As a consequence, research in nonlinear approximation theory has mostly focused on specific dictionaries. Mathematicians have studied specific dictionaries-spikes and sinusoids (Reference 4), wavelet packets (References 5, 25, and 24), frames (Reference 6), time/frequency tilings with algebraic hierarchical structure (Reference 19)- and presented individual algorithms (some provable, some heuristic) for constructing the best B-term approximation for any given function. Also, certain “mathematical recipes” such as algorithms for infinite dimensional Hilbert and Banach spaces have been proposed (see References 22 and 23). Save these exceptions, no provable, algorithmic results-running time for finding provably approximate B-term representations for inputs over general redundant dictionaries-are currently known.

[0039] 1.3 Local Search Heuristics

[0040] Most of the heuristic algorithms for general redundant dictionaries are local search methods. We briefly summarize them because they are relevant to our results. There are two common approaches.

[0041] 1. Find the j such that φ=φ_(j) maximizes (A−R, φ) over φ ε D. Update R by R←R+

A−R, φ

φ, and repeat a total of B times. More generally, find the set Λ of B′≦B of j's with the largest dot products, and put R←R+Σ_(jεΛ)

A, φ_(j)

φ_(j). We call this technique B′—fold matching pursuit, B′-MP, for B′≧1.

[0042] 2. Maintain a small subset Λ⊂D. Find the j such that φ=φ_(j) maximizes

A−R, φ

over φ ε D, and update Λ←Λ ∪ {φ_(j)}. Update R to be the optimal representation of the signal over the subspace spanned by {φ:φ ε Λ}, and repeat this process a total of B times. This technique is called orthogonal matching pursuit, OMP.

[0043] Both of these approaches have been studied extensively in the literature (MP appears in Reference 18 and OMP appears in Reference 21). However, there are two drawbacks.

[0044] The first issue is of proving bounds on the error of the solution. We note that in mathematical literature, the usual guarantee proved—if any bound at all is given-is that the error of the representation drops off at a prescribed rate as a function of B, for all signals in a given input class. From an algorithmic viewpoint such as ours, we wish to compare the error of the algorithms' output to the error of the optimal representation, whatever the optimal error may be, for a given signal and B. This viewpoint appears to be novel in mathematical approximation theory literature.

[0045] If the dictionary is an orthonormal basis, then the local search techniques above are equivalent and provably find the global optimum. In general, however, these methods do not provide any useful approximation results for finding the best representation for input functions (these are polynomial time heuristics for the NP-hard problem). Such an approximation result is not known even for special redundant dictionaries studied in the literature.

[0046] The second issue is of providing an efficient implementation, in particular, for performing each step of the local search. All known algorithms require at least Ω(|D|) time. Indeed, the general problem of finding the φ from the dictionary D that has the largest inner product with the signal is equivalent to the Farthest Neighbors Problem, which faces the same “dimensionality bottleneck” as the Nearest Neighbors Problem (see References 14, 15, and 17). Designing efficient algorithms for this problem is still an open challenge, in particular, since the local search heuristics rely on finding the optimal such φ's. This is a serious drawback because the community of researchers in highly nonlinear approximation theory are applied mathematicians who care about practical, implementable algorithms.

[0047] 1.4 Our Contributions

[0048] One might be tempted to study the problem assuming the dictionary is arbitrary so as to be general. However, an arbitrary dictionary is not a well-motivated choice. We would like to make the dictionaries as large as possible to achieve as high a rate of compression as possible; however, for specific classes of applications (e.g., encoding low-quality speech signals, compressing high motion video), the salient features of these signals are well-described by specific inherently suitable dictionaries. These dictionaries exhibit certain a naturalness for these applications such as smoothness, oscillatory behavior, etc. Redundant dictionaries must, therefore, balance the two needs-naturalness and succinctness of vectors on one hand and size on the other—and designing them is sometimes an art. We would like our algorithmic results to apply to largest number of redundant dictionaries proposed by domain experts in this area.

[0049] Further, it is natural to consider the scenario when the dictionary is provided ahead of time for preprocessing. The dictionary is typically large, say, of size Ω(N). Our goal is to process the input signal and determine a (near) optimal representation very efficiently. In particular, we would prefer to avoid scanning the dictionary entirely while processing the input signal.

[0050] We identify and focus on redundant dictionaries of small coherence. A μ-coherent dictionary D has coherence μ, for 0≦μ≦1, if |

φ₁, φ₂

|≦μ for all distinct φ₁, φ₂ ε D. A typical way to generate redundant dictionaries in practice is to take several orthonormal bases and combine them, and one gets the most out of this combination if the bases were as orthogonal to each other as possible (using sinusoids and spikes, as in Reference 11, is an example). Dictionaries generated in this manner have small coherence, thus motivating our work here. Coherence has been discussed as an important notion in the context of local search heuristics in the literature (see Reference 18), but we appear to the first to formalize the concept of coherent dictionaries and study the algorithmic problem of representing function near-optimally over μ-coherent dictionaries.

[0051] We make two contributions. Recall that the problem is, given input signal A and parameter B, determine the optimal B-term representation R over dictionary D of coherence μ, i.e., such that ∥A−R∥ is minimized.

[0052] We present an algorithm that for any ε, μB²=O(ε), finds an (1+ε)—factor approximation to the optimal representation in time N(B log(|D|)/ε)^(O(1)); the dictionary preprocessing takes time and space polynomial in |D|. We note that we assume that entries in A are bounded by N^(O(1)). The general result requires an additional factor in time and space polynomial in the number of bits in entries of A.

[0053] This is the first known provably approximate result for approximating a function under any redundant dictionary with small coherence O(ε/B²). In addition, it is very fast, taking time nearly linear in the input size modulo the polylogarithmic factor. This is an exponential speedup over previous heuristics that take time Ω( |D|) (recall that for redundant dictionaries, |D| dominates N). For the dictionary which consists of spikes and sinusoids, our approach specializes to give an algorithm that finds a nearoptimal representation in O(N²) time, improving the previous best (see Reference 4) of O(N^(3.5)) time via Linear Programming.

[0054] We explore the concept μ-coherent redundant dictionaries further, adding to the growing knowledge of redundant dictionary construction (see References 3, 10, and 11).

[0055] Known redundant dictionaries are of size only O(N) or O(N log(N)); it is desirable to construct larger dictionaries if “natural” properties of a specific application can still be captured. We consider image processing applications and propose a dictionary of significantly larger size: N ^(3/2)/B⁶, that we call segmentlets. This generalizes a number of known natural constructions in that area (see References 3, 10 and 11), such as beamlets, wedgelets, and ridgelets. By design, segmentlets are μ-coherent dictionaries for small μ and hence our main algorithmic result applies to them too.

[0056] Without focusing on any application, we focus on constructing large redundant dictionaries of small coherence. Using the Nisan-Wigderson (Reference 20) combinatorial design, we obtain such a dictionary of size exponential in N. We note that there are other combinatorial design that provide such large dictionaries with small coherence. However, this one can be constructed using small-space which may prove valuable for some application. Our motivation for including the Nisan-Wigderson combinatorial design was primarily to show that even though we focus on small conherence, nontrivial and exponentially large dictionaries can still be constructed.

[0057] The first category of result above is our main result, presented in §2. All our results on specific redundant dictionaries can be found in §3.

[0058] 1.5 Technical Approach.

[0059] Our algorithm is two-phased local search, OMP followed by MP. The crux of our technical contribution is our proof that in O(B/ε) iterations, our algorithm converges to within 1+ε of the optimal representation. At each iteration, this requires one to determine the dictionary element φ_(j) whose dot product has magnitude at least (1−η) times the largest dot product magnitude, and it suffices to estimate

A, φ_(j)

to within η∥A∥, additively, for some appropriate η. We are able to implement each such iteration using the approximate nearest neighbors data structures that have been designed recently. That gives the overall result.

[0060] There are some similarities between this work and our earlier work on finding near-optimal representations of functions over Fourier (see Reference 13) and Haar wavelet bases (see Reference 12), but also some crucial differences. Both Fourier and wavelet bases are non-redundant dictionaries, and so efficient algorithms already existed for finding even optimal representations over them. The emphasis in the previous work was on working under additional constraints: using only small (polylogarithmic) space and time, a 1a sampling or streaming. In our case here with redundant dictionaries, even given polynomial time and space, no previous result was known for optimal (or near-optimal) representation of functions. We do not focus on polylogarthmic space/time models, and leave that additional complication open. From a technical point of view, either OMP or MP by itself will suffice if we wanted to specialize our result for non-redundant basis such as Fourier or Haar wavelet; furthermore, one would not need approximate nearest neighbors. Thus, both the proof of approximation as well as the algorithms are more sophisticated for redundant dictionaries.

[0061] 2 Representation over Redundant Dictionaries

[0062] We consider vectors over R^(N) for simplicity; everything in this paper can be extended to C^(N) in an obvious way.

[0063] 2.1 Small Coherence Dictionary

[0064] DEFINITION 2.1. A set D={φ_(i)} of elements from R^(N) is a dictionary with coherence μ if span(D)=R^(N), ∥φ_(i)∥=1 for all i, and, for all distinct i and j, |

φ_(i), φ_(j)

|≦μ.

[0065] The definition of dictionary gives a condition only on pairs of dictionary elements. The following lemma gives an implication to larger sets of vectors.

[0066] LEMMA 2.1. Let φ₀, φ₁, . . . , φ_(B) be an arbitrary set of B+1 vectors from a dictionary of coherence μ. We have

[0067] 1. If μB<1, then {φ₁, . . . , φ_(B)} is independent.

[0068] 2. If μB<½, then the projection φ₀ onto the span of the other B vectors has norm at most {square root}{square root over (2 μ²B)}.

[0069] 3. If μB<½, then there exists a set {ψ_(i):i=1, . . . , B} of vectors such that:

[0070] The ψ's form an orthonormal system.

[0071] span(ψ₁, . . . , ψ_(B))=span{φ₁, . . . , φ_(B)}

[0072] ∥ψ_(i)−φ_(i)∥²≦8 μ²B.

[0073] Proof. First consider taking any linear combination φ₁, . . . , φ_(B). We have ∥Σ_(i)=1^(B)a_(i)φ_(i)∥²≧Σ_(i)=1^(B)a _(i) ²−Σ_(i≠j)|a_(i)a_(j)

φ_(i), φ_(j)

|, and Σ_(i≠j)|a_(i)a_(j)

φ_(i), φ_(j)

|≦μΣ_(i≠j)|a_(i)a_(j)|≦μ(Σ_(i)=1^(B)|a_(i)|)²≦μBΣ_(i)=1^(B)a_(i) ² by the Cauchy-Schwarz inequality, so $\begin{matrix} {{{\sum\limits_{i = 1}^{B}{a_{i}\varphi_{i}}}}^{2} \geq {\left( {1 - {\mu \quad B}} \right)\quad {\sum\limits_{i = 1}^{B}a_{i}^{2}}}} & (2.1) \end{matrix}$

[0074] This gives the first statement.

[0075] The length of the projection π₀ of φ₀ onto span(φ₁, . . . φ_(B)) is equal to the dot product

π₀,u₀

of π_(U) with the unit vector u_(U) along π_(U), which also equals

φ₀,u₀

. By the Cauchy-Schwarz inequality, this is the maximum, over all unit vectors u εspan(φ₁, . . . , φ_(B)), of

φ₀,u

. Write $u = \frac{\Sigma_{i}a_{i}\varphi_{i}}{{\Sigma_{i}a_{i}\varphi_{i}}}$

[0076] where Σ_(i)a_(i) ²=1; note that the denominator is at least {square root}{square root over (1−μB)}, by the above. Thus ∥π₀∥² is at most max _(∥u∥) _(²) ₌₁

u, φ₀

², which is at most ${{\max_{{\Sigma_{i}a_{i}^{2}} = 1}\frac{{\langle{{\Sigma_{i}a_{i}\varphi_{i}},\varphi_{0}}\rangle}^{2}}{{{\Sigma_{i}a_{i}\varphi_{i}}}^{2}}} \leq {\max_{{\Sigma_{i}a_{i}^{2}} = 1}\frac{{\langle{{\Sigma_{i}a_{i}\varphi_{i}},\varphi_{0}}\rangle}^{2}}{1 - {\mu \quad B}}} \leq {\max_{{\Sigma_{i}a_{i}^{2}} = 1}\frac{\left( {\Sigma_{i}\mu \quad a_{i}} \right)^{2}}{1 - {\mu \quad B}}} \leq \frac{\mu^{2}B}{1 - {\mu \quad B}} \leq {2\quad \mu^{2}B}},$

[0077] using the Cauchy-Schwarz inequality. This gives the second statement.

[0078] As for the third statement, let ψ′₁ be φ₁−π₁, where π₁ is the projection of φ₁ onto the span of {φ₂, . . . , φ_(B)}. Thus ψ′₁ is orthogonal to {φ₂, . . . , φ_(B)} and span(ψ′₁, φ₂, φ₃, . . . , φ_(B))=span(φ₁, φ₂, φ₃, . . . , φ_(B)).

[0079] By the above, ∥ψ′₁−φ₁∥²=∥π₁∥²≦2μ²B. Let $\psi_{1} = {\frac{\psi_{1}^{\prime}}{\psi_{1}^{\prime}}.}$

[0080] Since ∥ψ′₁∥²+∥π₁∥²=∥φ₁∥²=1, it follows that 1−2μ²B≦∥ψ′₁∥²≦1, so {square root}{square root over (1−2μ²B)}≦∥ψ′₁∥≦1, and ∥ψ₁−ψ′₁ ∥≦1−{square root}{square root over (1−2μ²B)}≦2μ²B. Since ∥ψ′₁−φ₁ ∥≦{square root}{square root over (2μ²B)}, it follows ∥ψ ₁−φ₁ ∥≦{square root}{square root over (2μ²B)}+2μ²B≦2{square root}{square root over (2μ²B)}, using the fact that x ²≦x for 0≦x≦1.

[0081] Recursively find ψ₂, . . . , ψ_(B) that are orthogonal, have the same span as φ₂, . . . , φ_(B), and such that ∥ψ_(i)−φ_(i)∥²≦8μ²(B−1)≦8μ²B.

[0082] Given as above {φ_(i)} ⊂ D, we say that {ψ_(i)} is an orthogonalization of {φ_(i)}.

[0083] 2.2 Overall Algorithm

[0084] Our overall algorithm is as follows. Starting with the zero representation R₁=0, we perform 0MP (the “error phase”) until ∥A−R₁∥²≦64B∥A−R_(opt)∥². Suppose that R₁ consists of B′<B terms. We then perform a single round of MP (the “energy phase”) to find a near-best (B−B′)—term representation R₂ to A−R₁. We then output R₁+R₂.

[0085] More quantitatively, we proceed as follows. Below, we first consider the energy phase. That is, given a signal A, we show how to find a representation R for A with square error ∥A−R∥²≦∥A−R_(opt)∥²+ε′∥A∥²; i.e., worse than optimal by a small multiple of the energy of the signal. Our ultimate goal is the stronger statement, ∥A−R∥²≦(1+ε)∥A−R_(opt ∥) ²=∥A−R_(opt)∥²+ε∥A−R_(opt)∥², where εand ε′ depend on B and μ. To do that, we show a modest error result (that is used first in the overall algorithm): given a signal A, we can find a representation R to A with square error ∥A−R∥²≦64B∥A−R_(opt)∥². That is, the square error of R is at most a moderate multiple of the optimal square error to the signal. Finally, combining these results, letting δ=∥A−R_(opt)∥ denote the optimal square error, we first reduce the square error to 64Bδ², then, by representing the residual, we reduce the square error to be additively suboptimal by ε′(64Bδ²). The result will have square error (1+64Bε′)δ²=(1+ε)δ², by definition of ε.

[0086] 2.3 Algorithm Details

[0087] 2.3.1 Energy Phase

[0088] We first show that we can roughly compare the error of two candidate representations for A, Σ_(iεΛ) ₁ a_(i)φ_(i)and by Σ_(iεΛ) ₂ β_(i)φ_(i) comparing Σ_(iεΛ) ₁

A, φ_(i)

with Σ_(iεΛ) ₂

A, φ_(i)

[0089] LEMMA 2.2. Suppose μB≦1. Let R₁=Σ_(iεΛ) ₁ a_(i)φ_(i) and R₂=Σ_(iεΛ) ₂ β_(i)φ_(i) be two B-term representations for A, such that R_(j) is the optimal representation in span (∪_(iεΛ) _(j) φ_(i)), j=1,2. If Σ_(iεΛ) ₁

A, φ_(i)

²≧Σ_(iεΛ) ₂

A, φ_(i)

², then

∥A−R ₁∥² ≦∥A−R ₂∥²+32μB∥A∥ ²

[0090] Proof. Let {ψ_(i):i ε Λ₁} be an orthogonalization of {φ_(i):i ε Λ₁}. Observe that R_(j)=Σ_(iεΛ) _(j)

A, ψ_(i)

ψ_(i), since each side is the unique best representation over span ({φ_(i):i ε Λ_(j)}=span({ψ_(i):i ε Λ_(j)}). Since ∥A−R_(j)∥²=∥A∥²−Σ_(iεΛ) _(j)

A, ψ_(i)

², it suffices to show that |Σ_(iεΛ) _(j)

A, ψ_(i)

²−Σ_(iεΛ) _(j)

A, φ_(i)

²|≦16μB∥A∥²

[0091] Proceeding, we have $\begin{matrix} {{\sum\limits_{i \in \bigwedge_{j}}{{A,\psi_{i}}}^{2}} = {\sum\limits_{i \in \bigwedge_{j}}{{{\langle{A,\varphi_{i}}\rangle} + {\langle{A,{\psi_{i} - \varphi_{i}}}\rangle}}}^{2}}} \\ {= {\sum\limits_{i \in \bigwedge_{j}}\left( {{{\langle{A,\varphi_{i}}\rangle}} - {{\langle{A,{\psi_{i} - \varphi_{i}}}\rangle}}} \right)^{2}}} \\ {= {{\sum\limits_{i \in \bigwedge_{j}}{{\langle{A,\varphi_{i}}\rangle}}^{2}} - {2\quad {\sum\limits_{i \in \bigwedge_{j}}{{{\langle{A,\varphi_{i}}\rangle}}{{{\langle{A,{\psi_{i} - \varphi_{i}}}\rangle}}.}}}}}} \end{matrix}$

[0092] A bound for the last term will be reused below, so we isolate it. By Lemma 2.1, it is at most 2{square root}{square root over (8μ²B)}∥A∥Σ_(iεΛ) _(j) |

A, φ_(i)

| which, by the Cauchy-Schwarz inequality, is at most 2{square root}{square root over (8μ²B)}∥A∥{square root}{square root over (BΣ_(iεΛ) _(j) |

A, φ_(i)

|²)}. Continuing, $\begin{matrix} {\sqrt{B\quad {\sum\limits_{i \in \bigwedge_{j}}{{\langle{A,\varphi_{i}}\rangle}}^{2}}} \leq \sqrt{B{\sum\limits_{i \in \bigwedge_{j}}\left( {{{\langle{A,\psi_{i}}\rangle}} + {{\langle{A,{\varphi_{i} - \psi_{i}}}\rangle}}} \right)^{2}}}} \\ {\leq \sqrt{2B\quad {\sum\limits_{i \in \bigwedge_{j}}\left( {{{\langle{A,\psi_{i}}\rangle}}^{2} + {{\langle{A,{\varphi_{i} - \psi_{i}}}\rangle}}^{2}} \right)}}} \\ {{\leq \sqrt{{2B{A}^{2}} + {8\quad \mu^{2}B^{3}\quad {A}^{2}}}},} \end{matrix}$

[0093] so that 2Σ_(iεΛ) _(j) |

A, φ_(i)

∥

A, ψ_(i)−φ_(i)

|≦8μB∥A∥² {square root}{square root over (1+4μ²B²)}≦8μB{square root}{square root over (2)}∥A∥², and Σ_(iεΛ) _(j)

A, ψ_(i)

²≧Σ_(iεΛ) _(j)

A, φ_(i)

²−16μB∥A∥²

[0094] Similarly, Σ_(iεΛ) _(j) |

A, ψ_(i)

|²≦Σ_(iεΛ) _(j) |

A, φ_(i)

|²+2Σ_(iεΛ) _(j) |

A, φ_(i)

|·|

A, ψ_(i)−φ_(i)

|+Σ_(iεΛ) _(j) |

A, ψ_(i)−φ_(i)

|², which is at most Σ_(iεΛ) _(j) |

A, φ_(i)

|²+8μB{square root}{square root over (2)}∥A∥²+8μ²B²∥A∥²≦Σ_(iεΛ) _(j) |

A, φ_(i)

|²+16μB∥A∥², since μB≦{fraction (1/2)}. The result follows.

[0095] Algorithmically, one finds the B vectors ø ε D with largest dot products to A, then finds the best representation R to A over the span of those vectors (for example, by orthogonalizing the vectors first). The resulting representation R will have square error ∥A−R∥²≦∥A−R_(opt)∥²+ε∥A∥², for ε=32μB.

[0096] COROLLARY 2.1. For a signal A, a single iteration of B-fold Matching Pursuit over a incoherent dictionary D returns a representation R with ∥A−R∥²≦∥A−R_(opt)∥²+32μB∥A∥²

[0097] 2.3.2 Error Phase

[0098] Lemma 2.4 below says that if we have a representation R whose error is significantly worse than optimal (roughly, error {square root}{square root over (Bδ)} compared with optimal value of δ), then the ø with largest dot product represents so much of the signal that we are forced to take it in any optimal or near-optimal representation. This principle will be used, later, to show correctness of a greedy algorithm.

[0099] First we prove a lemma that will be useful.

[0100] LEMMA 2.3. Fix B>0. Let Λ be a signal. Let A be a set of size less than B, and let R=Σ_(iεΛ)a_(i)φ_(i) be the best representation for A over {ø_(i):i ε Λ}. Let R_(opt)=Σ_(iεΛ) _(opt) a_(i) ^(opt)φ_(i) be the best B-term representation for A over all D), subject to Λ⊂Λ_(opt). (Note that ai does not necessarily equal a_(i) ^(opt).) If ∥A−R∥²≦(1+ε)∥A−R_(opt)∥² for 64μ²B² ≦ε<1, then there exists i ε Λ _(opt)/Λsuch that ${\langle{{A - R},\varphi_{i}}\rangle}^{2} \geq {\frac{ɛ}{8\quad B}{{A - R}}^{2}}$

[0101] Proof. Note that (1+ε)>(1−ε/2)⁻¹, so the hypothesis ∥A−R∥²≧(1+ε)∥A−R_(opt)∥² implies (1−ε/2)∥A−R∥²≧∥A−R_(opt)∥², so that $\begin{matrix} {{{{A - R}}^{2} - {\frac{ɛ}{2}{{A - R}}^{2}}} \geq {{{A - R_{opt}}}^{2}\quad {and}}} \\ {{{{A - R}}^{2} - {{A - R_{opt}}}^{2}} \geq {\frac{ɛ}{2}{{{A - R}}^{2}.}}} \end{matrix}$

[0102] Let {ψ_(i)} be an orthogonalization of {φ_(i):φ_(i) ε Λ_(opt)} extending an orthogonalization of {φ_(i):φ_(i) ε Λ} Then R_(opt) is the best representation over {ψ_(i):i εΛ_(opt)} and R is the best representation over {ψ_(i):i ε Λ}. Furthermore, by orthogonality, R_(opt)=Σ_(iεΛ) _(opt)

A, ψ_(i)

ψ_(i) and R=Σ_(iεΛ)

A, ψ_(i)

ψ_(i), using the same coefficients as R_(opt).

[0103] By Parseval's equality and orthogonality,

∥A−R∥ ² −∥A−R _(opt)∥² =∥R _(opt−) R∥ ²

[0104] which equals Σ_(iεΛ) _(opt) _(\Λ)

A, ψ_(i)

², So it follows that, for some i ε Λ_(opt)\Λ, we have ${\langle{A,\psi_{i}}\rangle}^{2} \geq \frac{{{A - R}}^{2} - {{A - R_{opt}}}^{2}}{B} \geq {\frac{ɛ}{2B}{{{A - R}}^{2}.}}$

[0105] Since R is orthogonal to ψ_(i), it follows that ${\langle{{A - R},\psi_{i}}\rangle}^{2} = {{\langle{A,\psi_{i}}\rangle}^{2} \geq {\frac{ɛ}{2B}{{A - R}}^{2}}}$

[0106] Now, ${{{\langle{{\varphi_{i} - \psi_{i}},{A - R}}\rangle}} \leq {{{\varphi_{i} - \psi_{i}}} \cdot {{A - R}}} \leq {\left( \sqrt{8µ^{2}B} \right){{A - R}}} \leq {\left( \sqrt{\frac{ɛ}{8B}} \right){{A - R}}}},$

[0107] so ${{{\langle{{A - R},\varphi_{i}}\rangle}} \geq {{{\langle{{A - R},\psi_{i}}\rangle}} - {{\langle{{A - R},{\psi_{i} - \varphi_{i}}}\rangle}}} \geq {\left( {\sqrt{\frac{ɛ}{2B}} - \sqrt{\frac{ɛ}{8B}}} \right){{A - R}}}} = {\sqrt{\frac{ɛ}{8B}}{{A - R}}}$

[0108] , and the result follows.

[0109] LEMMA 2.4. Let A be a signal and let R=Σ_(iεΛ)α_(i)φ_(i) be the best representation for A over Λ of size less than B. Suppose there's a set A_(opt) ⊃Λ with |A_(opt)|=B, and a representation R_(opt)=Σ_(iεΛ) _(opt) β_(i)φ_(i) for A such that ∥A−R∥²>64B∥A−R_(opt)∥². Finally, assume ${8{µ^{2}\left( {B + 1} \right)}} < {\frac{1}{64\quad B}.}$

[0110] If ø=ø_(i) maximizes (A−R, ø) over ø ε D, then i ε Λ_(opt.)

[0111] Proof. Suppose not. By Lemma 2.3 with ε=½, since φ_(i) is maximal, ${\langle{{A - R},\varphi_{i}}\rangle}^{2} \geq {\frac{1}{16\quad B}{{{A - R}}^{2}.}}$

[0112] Let {ψ_(j)} denote an orthogonalization of {φ_(j):j ε Λ_(opt) ∪ {i}}.

[0113] Then $\begin{matrix} {{\langle{{A - R},\psi_{i}}\rangle} \geq {{\langle{{A - R},\varphi_{i}}\rangle} - {{{A - R}}{{\varphi_{i} = \psi_{i}}}}} \geq} \\ {{{\left( {\sqrt{\frac{1}{16\quad B}} - \sqrt{8{µ^{2}\left( {B + 1} \right)}}} \right){{A - R}}} \geq}} \\ {{{\left( {\sqrt{\frac{1}{16\quad B}} - \sqrt{\frac{1}{64\quad B}}} \right){{A - R}}},}} \end{matrix}$

[0114] so that ${\langle{{A - R},\psi_{i}}\rangle}^{2} \geq {\frac{1}{64\quad B}{{A - R}}^{2}} > {{{A - R_{opt}}}^{2}.}$

[0115] On the other hand, since ψ_(i) is orthogonal to Λ_(opt) ∃ span(Λ_(opt)\Λ), it follows that

A−R, ψ_(i)

=

A−R_(opt), ψ_(i)

+

R_(opt)−R, ψ_(i)

=

A−R_(opt), ψ_(i)

Also, by optimality of R_(opt), A−R_(opt) is orthogonal to R_(opt), so that ∥A−R_(opt)∥²=∥A∥²−∥R_(opt) ∥. Finally, since ψ_(i) is orthogonal to Λ_(opt) ⊃Λ, Λ∪{i} can be extended to a basis, and we have ∥A∥²≧∥R_(opt)∥²+

A−R_(opt), ψ_(i)

², i.e.

A−R, ψ_(i)

₂=

A−R_(opt), ψ_(i)

²≦∥A∥²−∥R_(opt)∥²=∥A−R_(opt)∥², a contradiction.

[0116] Algorithmically, we can start with a signal A and representation R=0 over subspace Λ=Ø. As long as ∥A−R∥²>64B∥A−R_(opt)∥², we can identify a vector φ εD that is in every optimal representation. We add φ to Λ, let R be the best representation over A, and continue looping. We terminate after at most B iterations, and, when we terminate, ∥A−R∥≦64B∥A−R_(opt)∥.

[0117] COROLLARY 2.2. For a signal A, Orthogonal Matching Pursuit, in at most B iterations, over a p-coherent dictionary D, returns a representation R with ∥A−R∥²≦64B∥A−R_(opt)∥². OMP stopped when ∥A−R∥²<64B∥A−R_(opt)∥² returns a representation over a subspace of dimension B′≦B that can be extended to a B-dimensional subspace containing an optimal representation.

[0118] 2.3.3 Putting it All Together

[0119] THEOREM 2.1. Fix a dictionary D with coherence μ. Let A be a signal and suppose it has a B-term representation over D with error ∥A−R_(opt)∥=δ, where B<1/(32μ). Then, in iterations polynomial in B, we can find a representation with error at most {square root}{square root over ((1+2064μB²)δ)}.

[0120] Proof. The algorithm is as follows. Assume we know δ=∥A−R_(opt)∥ in advance; we will remove this assumption below. Use Corollary 2.2 to find a B′-term representation R₁ over Λ⊂D with ∥A−R∥²≦64Bδ² and such that Λ₁ is a subset of a space containing an optimal representation. Then use Corollary 2.1 to find a representation R₂·with square error at most 32μB∥A−R∥²≦32·64μB²∥A−R_(opt)∥² worse than the best representation for A−R₁. Output R₁+R₂.

[0121] We note that the hypotheses on μ and B in Corollaries 2.2 and 2.1 are implied by our hypothesis B<1/(32μ). Since, as we show below, the best representation for A−R₁ has square error at most δ²(1+512μ²B³), the overall representation has square error at most δ²(1+512μ²B³+2048μB²)≦(1+2064μB²)δ², since μB<{fraction (1/32)}. It remains only to show that A−R₁ has a representation with square error at most δ²(1+512μ²B³).

[0122] To see this, fix an optimal representation R_(opt)=Σ_(iεΛ)β_(i)φ_(i) consisting of vectors in some Λ⊃Λ₁. Let {ψ_(i)} denote an orthogonalization of {φ_(i):φ_(i)εΛ} that extends an orthogonalization of {φ_(i):φ_(i)εΛ₁}. Then R_(opt)=Σ_(iεΛ)

A, ψ_(i)

ψ_(i) and R₁=Σ_(iεΛ) ₁

A, ψ_(i)

ψ_(i).

[0123] Consider R_(opt)−R₁, which is orthogonal to R₁ and to A−R_(opt). We first claim that R_(opt)−R has some good (B−|Λ₁|)—term representation over D. Specifically, we will approximate R_(opt)−R₁=Σ_(iεΛ\Λ) ₁

R_(opt)−R₁ψ_(i)

ψ_(i) by Σ_(iεΛ\Λ) ₁

R_(opt)−R₁ψ_(i)

φ_(i). Since each ∥ψ₁−φ₁∥ is small, we expect that substituting φ_(i) for ψ_(i) in the expansion for R_(opt)−R to be a small perturbation.

[0124] We have ∥(R_(opt)−R₁)−Σ_(iεΛ\Λ) ₁

R_(opt)−R₁, ψ_(i)

φ_(i)∥²=∥Σ_(iεΛ\Λ) ₁

R_(opt)−R₁, ψ_(i)

(ψ_(i)−φ_(i))∥² is at most Σ_(iεΛ\Λ) ₁

R_(opt)−R₁, ψ_(i)

²Σ_(iεΛ\Λ) ₁ ∥ψ_(i)−φ_(i)∥², which is at most $\begin{matrix} {{{{{R_{opt} - R_{1}}}^{2} \cdot {B\left( {8µ^{2}B} \right)}} \leq}} \\ {{{{{A - R_{1}}}^{2} \cdot {B\left( {8µ^{2}B} \right)}} \leq}} \\ {{{64\quad B\quad {\delta^{2} \cdot {B\left( {8µ^{2}B} \right)}}} = {512\quad µ^{2}B^{3}{\delta^{2}.}}}} \end{matrix}$

[0125] Since A−R_(opt) is orthogonal to span(A), it follows that the Σ_(iεΛ\Λ) ₁

R_(opt)−R₁ψ_(i)

φ_(i) gives a representation for A−R₁ with corresponding error, namely $\begin{matrix} {{{\left( {A - R_{1}} \right) - {\sum\limits_{i \in {\Lambda \smallsetminus \Lambda_{1}}}{{\langle{{R_{opt} - R_{1}},\psi_{i}}\rangle}\varphi_{i}}}}}^{2} = {{{A - R_{opt}}}^{2} +}} \\ {{{\left( {R_{opt} - R_{1}} \right) - {\sum\limits_{i \in {\Lambda \smallsetminus \Lambda_{1}}}{{\langle{{R_{opt} - R_{1}},\psi_{i}}\rangle}\varphi_{i}}}}}^{2}} \\ {= {\delta^{2}\left( {1 + {512µ^{2}B^{3}}} \right)}} \end{matrix}$

[0126] Using Lemma 2.2, one can find a (B−|Λ₁|)-term approximation R₂ to A−R₁ whose error is at most 32μB∥A−R₁∥²≦32μB(64Bδ²) worse, additively, then the best such representation. It follows that A − (R₁ + R₂)² ≤ (A − R₁) − R₂² ≤ δ²(1 + 512  µ²B³) + 32  µ  B(64  B  δ²) ≤ δ²(1 + 16  µ  B² + 2048  µ  B²),  since  µ  B < 1/32 ≤ δ²(1 + 2064  µ  B²)

[0127] Now suppose we do not know δ in advance. Then we would not know B′=|Λ₁|, i.e., we would not know when to switch from the error phase to the energy phase. We simply try all possibilities for B′<B.

[0128] So far, we have assumed that one can find the ø ε D with largest |

A, ø

| and estimate that dot product in unit time. We can, in fact, make this algorithm feasible and obtain a result similar to Theorem 2.1 using Nearest Neighbor Data Structures (see References 14, 15, and 17). For some c, fix ε, cμB²≦ε≦1; our goal will be to find a representation with error at most (1+ε) worse than optimal.

[0129] For Lemmas 2.2-2.4 and Theorem 2.1, we can prove that it suffices to find ø with near-biggest dot product and to approximate the dot product; i.e., to find ø_(j) such that, for all i, |

A, φ_(j)

|²≧|

A, φ_(i)

|²−η∥A∥² and to estimate

A, φ_(j)

as

A, φ_(j)

with |

A, φ_(j)

−

A, φ_(j)

|²≦η∥A∥², for some appropriate η, polynomially related to B/ε. There is a limited set S of vectors for which we will want dot products; specifically, we will want (ψ₁ψ₂) only if each ψ is a linear combination of at most B vectors from D ∪ {A}, in which the coefficients are, by (2.1), at most O(∥A∥²) (here taken to be N^(O(1))) and written to the unit of precision (here taken to be 1). Thus there are at most (|D|+1)^(O(B)) possible vectors in S. Therefore, we can use the following steps:

[0130] Normalize the signal to $\frac{A}{A},$

[0131]  find a representation for the normalized signal, and scale back up. For the normalized signal, we can get the approximate dot product operations from approximate l² distance operations.

[0132] Using Reference 1, randomly project the signal A and the dictionary vectors, using a randomly-chosen linear function h from a particular family, into a log |S|/η^(O(1)) dimensional space so the l² norm between any pairs of vectors from S is approximated to factor (1+η).

[0133] Use the approximate nearest neighbor results in References 14, 15 and 17 to return a (1+η)—approximation to the closest dictionary vector in l² norm to a query of the form A−R, in time (Blog|D|/η)^(O(1)). This allows us to compute RI of Theorem 2.1. To compute R₂ of Theorem 2.1, we need to find the B largest dot products with A−R₁. To do that, repeatedly find the ø_(i) with largest dot product to A−R and, by properties of the nearest neighbor data structure, remove h(ø_(i)) from the dictionary, in time (Blog|D|/η)^(O(1)). Finally, to approximate the best representation of A over Λ, find the best formal representation of h(A) over {h(ø):ø ε Λ} and use those coefficients.

[0134] Note that, even for ø_(j) ε Λ, it is possible that

A−R, φ_(j)

·0. In fact,

A−R, φ_(j)

may be the maximum dot product over ø ε D, so we may need to choose oj on multiple iterations. One can show, however, that only a small number of iterations is needed for OMP to find a B′-term representation RI, B′<B, with ∥A−R₁∥²≦64B∥A−R_(opt)∥². The energy phase takes just a single iteration of (B−B′)-fold MP. In summary,

[0135] THEOREM 2.2. Fix a dictionary, D, of coherence A, over vectors of length N. For certain constants c and c′, fix B<c/μ and fix e, c′μB²<ε<1. There is an algorithm that preprocesses D taking time and space (B|D|/ε)^(O(1)). For any given signal, the algorithm produces a representation for A with error (1+ε) times that of the optimal representation taking time (B log|D|/ε)^(O(1)) only.

[0136] Note that some cost dependence on N is hidden by dependence on |D|≧N.

[0137] The idea to use nearest neighbor data structures in this context was suggested independently by Piotr Indyk.

[0138] 3. Special Dictionaries

[0139] In this section, we briefly explore the notion of small coherence dictionaries further. We consider several specific dictionaries and analyze their properties.

[0140] Three desired properties of redundant dictionaries are their large size, their naturalness, and the possibility to find representations over them quickly. Here “naturalness” encompasses the idea that signals typical of a particular application should have concise, low-error representations, and that the dictionary elements themselves should have semantic significance within the application. These goals are somewhat in conflict with each other; nevertheless, we show that our criterion of low coherence applies to several large and natural dictionaries in the literature (or small variations).

[0141] 3.1 Spikes and Sinusoids

[0142] Two fundamental orthonormal bases are the spikes and sinusoids. A spike is a vector δ_(s) defined by δ_(s)(s)=1 and δ_(s)(t)=0 for t≠s. A sinusoid is a complex-valued vector ψ_(w) defined by ${\psi_{w}(t)} = {\frac{1}{\sqrt{N}}{^{2\quad \pi \quad \quad {{wt}/N}}.}}$

[0143] It is immediate that the dot product of any spike with any sinusoid equals 1/{square root}{square root over (N)} in absolute value and the dictionary D formed by their union is μ-coherent, for μ=1/{square root}{square root over (N)}. It follows that the algorithm in Theorem 2.1 applies to this dictionary. Note that, in particular, if the signal is exactly represented as a sum of B spikes and sinusoids, so that the optimal error is zero, then the error of our output is zero-we recover the signal exactly.

[0144] We note that Theorem 2.1 gives a running time of O(B⁴+B²N+BN log(N)) to find a representation with error at most (1+O(B²/{square root}{square root over (N)})) times optimal. For each of B iterations, we have a subset Λ⊂D of spikes and sinusoids available for the representation. To find the best representation R over Λ, we could find an orthonormal basis for span(Λ), which can be done in time B³. We then need to find the ø_(j)'s with biggest dot products to A−R. We can find all

φ_(j), A

's explicitly in time O(N log(N)) and all

φ_(j), R

's symbolically in time BN. For the extreme case of B=O({square root}{square root over (N)}), our algorithm takes time O(N²) to recover exactly any signal that is exactly represented as B spikes and sinusoids. By contrast, in Reference 4, the authors give an algorithm with runtime O(N^(3.5)) which recovers the signal exactly, assuming it consists of B spikes and/or sinusoids. Even for the limited case of an exact representation, their algorithm is quite expensive.

[0145] Our algorithm and Reference 4 apply to the combination of any two incoherent orthonormal bases. The analysis above for our algorithm assumes that one of the bases is presented in the coordinate system defined by the other (equivalently, that one can take the dot product between any two dictionary elements in constant time.) If this is not the case, then one could compute and store all O(N²) such dot products in the time to multiply two N×N matrices, i.e., N^(2+a) for some a, 0<a<1.

[0146] Other speedups are possible for our algorithm, using Theorem 2.2 in general or using techniques in Reference 13 for the spikes and sinusoids dictionary in particular. Details will be given in the final version of this paper.

[0147] 3.2 Combinatorial Designs

[0148] The dictionary in the previous subsection had small coherence and was a combination of two fundamental orthonormal bases but it was not very large. In fact, we can build a dictionary with small coherence that is considerably larger using a combinatorial design, such as the matrix at the center of the Nisan-Wigderson generator. Specifically, in Reference 20, the authors show how to build a collection S of subsets of Reference 1, N, for N an even power of a prime, such that each subset has size exactly {square root}{square root over (N)}, any two subsets intersect in at most d places, and the number of subsets is Ω(N^((d+1)/2)).

[0149] Define D by D={N^(1/4)χs:SεS}. It follows that |D|=|S|=Ω(N^((d+1)/2)) and, for distinct φ₁ and φ₂ in D, |

φ₁, φ₂

|≦d1{square root}{square root over (N)}.

[0150] For d=2, we get a superlinear-sized dictionary of size roughly N^(3/2) and coherence roughly 2/{square root}{square root over (N)}. At the other extreme, one can take d almost as large as {square root}{square root over (N)}, yielding coherence less than 1 and |D| roughly 2^({square root}{square root over (N)}). One can also use values for d between the two extremes—for example, if d=log(N) or d=N^(1/4), the dictionary has N^(Ω(log(N)))>N^(O1) or 2^(Ω(N) ^(1/4)) >>N^(O(1)) elements and coherence log (N)/{square root}{square root over (N)} or N^(1/4), not much worse than the expected absolute dot product 1/{square root}{square root over (N)} of a pair of random unit vectors.

[0151] This construction gives us non-trivial dictionaries with extremely small coherence, extremely large size, or both, to which we can apply Theorem 2.1 and obtain efficient algorithms.

[0152] 3.3 Segmentlets

[0153] Next, we present a redundant dictionary that is inspired by beamlets (see Reference 11), a dictionary used for edge detection in image analysis. Consider the space of functions on the square array of side {square root}{square root over (N)}. Fix a parameter p, to be determined later. Consider the set of all line segments with endpoints (x₁, y₁) and (X₂, y₂), such that x₁, y₁, and Y₂ are multiples of p and X₂=x₁+p. Then the number of segments is Ω(N^(3/2)/p³), much greater than the dimension, N, of the space of functions, and any two segments intersect at most once and have horizontal extent exactly p.

[0154] From each line segment, we next construct a set of pixels, each of which intersects the line segment. In general, a line segment may intersect several contiguous pixels in the same column; in that case, put only the middle one or two intersected pixels into the associated set. It follows that any two sets of pixels constructed this way will intersect in at most O(1) pixels and each will consist of O(p) pixels. For the dictionary, take normalized characteristic functions on theses sets of pixels. The coherence is O(1/{square root}{square root over (p)}) and the size of the dictionary equals the number of segments, N^(3/2)/p³.

[0155] Suppose we are interested in B-term representations. To apply the above techniques, we need coherence less than 1/B, i.e., p≧B². It follows that the dictionary can have size N^(3/2)/B⁶, which is greater than N for sufficiently small B. In particular, the size of the dictionary is superlinear, by more than log factors. Thus the size of the segmentlet dictionary is larger than the beamlet dictionary, and we can apply Theorem 2.1 to obtain efficient algorithms for near-optimal B-term representations over segmentlets. Segmentlets are natural for capturing edges in images. We believe that this dictionary will have exciting applications in image processing.

[0156] 4. Remarks

[0157] We have presented algorithmic results for the basic problem in highly nonlinear approximation: An efficient algorithm for near-optimal representation of input signals over μ-coherent dictionaries. There are several additional insights we can provide into highly nonlinear approximation problems. We briefly describe them here.

[0158] The problem of representing functions using redundant dictionaries has a few variants. As defined above; a representation R for input A ε R^(N) is a linear combination of dictionary elements, R=Σ_(iεΛ)a_(i)φ_(i) for ø_(i) ε D. In general, there are two measures to assess the goodness of the representation. The first is how well it approximates A and we measure this error ∥A−R∥_(p) in l^(P) norm for 1≦p≦∞; the results in this paper focus on p=2. The second metric attempts to capture varying notions of parsimony and we measure this with the tq norm of the representation itself, ∥R∥_(q)=Σ_(iεΛ)|a_(i)|^(q) for q=0 or 1≦q≦∞. We focused throughout on minimizing the l² norm of the error subject to the l⁰ norm of the representation equaling zero. In general, other (p,q) combinations may be of interest, and one may want to fix the l^(P) norm of the error (for example to zero) and minimize the l^(q) norm of the representation. For example, instead of seeking the best B-term representation, we may seek the representation minimizing Σ_(i)|a_(i)|, i.e. ∥a∥₁. We have a polynomial time solution for this case via linear programming when p=∞. Let Φ be the matrix of vector ø_(i) at components j, Φ_(ij)=φ_(i)(j). Express each coefficient a_(i) as a sum of positive and negative parts, a_(i)=p_(i)−n_(i). The linear program is to minimize Σ_(i)(p_(i)+n_(i)) subject to $\left\{ \begin{matrix} {{\left( {\Phi - \Phi} \right)\begin{pmatrix} p \\ n \end{pmatrix}} \leq {A + ɛ}} \\ {{\left( {\Phi - \Phi} \right)\begin{pmatrix} p \\ n \end{pmatrix}} \geq {A - ɛ}} \end{matrix}\quad \right.$

[0159] Another variation is one in which we wish to minimize the energy of the representation, i.e., minimize ∥a∥₂. This problem has a polynomial time solution via semidefinite programming. This is because

minimize a ^(t) a subject to ∥Φa−A∥ _(p)≦ε

[0160] is a semidefinite program. Note that a^(t)a is the squared l² norm of the representation and the constraints define a convex region about the point A ε R^(N). We have initiated the study of the algorithmic complexity of highly nonlinear approximation problems. Approximation theory is a rich area, and it is active in mathematics, signal processing and statistics; theoretical computer scientists can have a significant impact here. Many algorithmic problems remain open. Finally, as remarked earlier, applied mathematicians in approximation theory care about practical algorithms. We believe that our algorithms here are not difficult to implement, but we leave that study for the future.

[0161]FIG. 1 is a block diagram of an exemplary embodiment of a telecommunications system 1000 that can implement an exemplary embodiment of a method and/or algorithm described herein. System 1000 can comprise any number of signal sources 1100. Any signal source 1110, 1120, 1130 can be an information device. From any signal source 1110, 1120, 1130 can emanate a signal 1112, 1122, 1132, respectively.

[0162] Any of the signal sources 1100 can be coupled to a network 1200. Coupled to network 1200 can be any number of information devices 1300 that can receive a signal and/or process a signal. The signal can be an analog or a digital signal. The content of the signal can be, for example, data, audio, speech, and/or video, etc. The signal can be carried by and/or modulated upon a wireless and/or wireline media, including an optical media. For example, the signal can be an optical voice communication; a radar and/or sonar signal; a microwave data transmission; an FM radio audio signal; a landline, cellular, or IP telephone speech signal; a satellite television signal; an e-mail or HTTP message sent over a DSL connection; and/or a Bluetooth data signal, etc.

[0163] A signal receiver 1310 can provide a received signal to a signal processor 1320 that can iteratively determine a provable near-optimal B-term representation of the signal. The determined representation can be useful for a wide variety of functions, such as encoding, modulating, compressing, encrypting, detecting, decrypting, decompressing, demodulating, and/or decoding, etc., the signal. Each B-term representation can be formed from a linear combination of a plurality of elements selected from a redundant dictionary. The dictionary can be predetermined and/or μ-coherent. The determined near-optimal B-term representation can correspond to an error that is a minimum of a set of errors corresponding to all possible B-term representations associated with the predetermined redundant μ-coherent dictionary.

[0164] Coupled to information device(s) 1300 can be a memory device 1400 that can store one or more dictionaries and/or possible B-term representations. Coupled to network 1200 can be an information device 1500 that can receive the near-optimal B-term representation and/or the B-term representation having the minimum error. Information device 1500 can render any received B-term representation on an output device 1510. Via a user interface 1520 of information device 1500 inputs can be received and/or outputs rendered regarding the dictionaries, errors, representations, and/or the values for B, μ, and/or other relevant parameters.

[0165]FIG. 2 is a flow diagram of an exemplary embodiment of a method 2000 for automatically implementing at least a portion of the nonlinear approximation algorithm described herein. At activity 2100, a signal is received. At activity 2200, a value for B is determined, obtained, and/or received. At activity 2300, a dictionary is determined, obtained, and/or received. At activity 2400, a B-term representation of the signal is determined. At activity 2500, an error for the determined B-term representation is determined. At activity 2600, if the determined error is a minimum compared to the errors and/or likely errors of the other possible B-term representations, then flow proceeds to activity 2700; otherwise, flow returns to activity 2400. At activity 2700, the B-term representation associated with the minimum error, sometimes referred to herein as the “near-optimal B-term representation”, is output and/or stored.

[0166] In certain exemplary embodiments, the near-optimal B-term representation is determined in a time proportional to N raised to a power of approximately 2, where N is a total number of unit vectors comprised by the predetermined redundant μ-coherent dictionary. In certain exemplary embodiments, the near-optimal B-term representation is determined in O(N²) time, where N is a total number of unit vectors comprised by the predetermined redundant μ-coherent dictionary.

[0167] In certain exemplary embodiments, the near-optimal B-term representation determined without entirely searching the predetermined redundant μ-coherent dictionary.

[0168] In certain exemplary embodiments, B is less than 1/(32μ). In certain exemplary embodiments, in O(B/ε) iterations an associated intermediate representation is within an error, denoted I+ε, of the near-optimal B-term representation.

[0169] In certain exemplary embodiments, any portion of method 2000 can be repeated in any defined manner, including periodically, pseudo-randomly, and randomly. In certain exemplary embodiments, any portion of method 2000 can occur dynamically.

[0170]FIG. 3 is a block diagram of an exemplary embodiment of an information device 3000, which in certain operative embodiments can represent, for example, signal source 1100, information device 1300, and/or information device 1500 of FIG. 1. Information device 3000 can comprise any of numerous well-known components, such as for example, one or more network interfaces 3100, one or more processors 3200, one or more memories 3300 containing instructions 3400, one or more input/output (I/O) devices 3500, and/or one or more user interfaces 3600 coupled to I/O device 3500, etc.

[0171] In certain exemplary embodiments, via one or more user interfaces 3600, such as a graphical user interface, a user can implement, facilitate, monitor, and/or manage an exemplary embodiment of the nonlinear approximation algorithm.

[0172] Still other embodiments will become readily apparent to those skilled in this art from reading the above-recited detailed description and drawings of certain exemplary embodiments. It should be understood that numerous variations, modifications, and additional embodiments are possible, and accordingly, all such variations, modifications, and embodiments are to be regarded as being within the spirit and scope of the appended claims. For example, regardless of the content of any portion (e.g., title, field, background, summary, abstract, drawing figure, etc.) of this application, unless clearly specified to the contrary, there is no requirement for the inclusion in any claim of the application of any particular described or illustrated activity or element, any particular sequence of such activities, or any particular interrelationship of such elements. Moreover, any activity can be repeated, any activity can be performed by multiple entities, and/or any element can be duplicated. Further, any activity or element can be excluded, the sequence of activities can vary, and/or the interrelationship of elements can vary. Accordingly, the descriptions and drawings are to be regarded as illustrative in nature, and not as restrictive. Moreover, when any number or range is described herein, unless clearly stated otherwise, that number or range is approximate. When any range is described herein, unless clearly stated otherwise, that range includes all values therein and all subranges therein. Any information in any material (e.g., a United States patent, United States patent application, book, article, etc.) that has been incorporated by reference herein, is only incorporated by reference to the extent that no conflict exists between such information and the other statements and drawings set forth herein. In the event of such conflict, including a conflict that would render a claim invalid, then any such conflicting information in such incorporated by reference material is specifically not incorporated by reference herein.

REFERENCES

[0173] The following references are incorporated by reference herein in their entirety:

[0174] 1. N. Alon, Y. Matias, and M. Szegedy. The space complexity of approximating the frequency moments. JCSS, 58(1):137-147, 1999.

[0175] 2. Constructive Approximation. On the world wide web at math.vanderbilt.edu/˜ca.

[0176] 3. E. Candes. Ridgelets: Theory and Applications. PhD thesis: Dept. of Statistics, Stanford University, 1998.

[0177] 4. S. S. B. Chen, D. L. Donoho, and M. A. Saunders. Atomic decomposition by basis pursuit. SIAM J. Sci. Computing, 20:33-61, 1999.

[0178] 5. R. Coifman and M. V. Wickerhanser. Entropy-based algorithms for best basis selection. IEEE 7Fans. Inform. Theory, 38(2), March 1992.

[0179] 6. I. Daubechies. Ten lectures on wavelets. SIAM, Philadelphia, 1992.

[0180] 7. G. Davis_S. Mallat and M. Avellaneda. Greedy adaptive approximation. Journal of Constructive Approximation, 13:57-98, 1997.

[0181] 8. R. A. DeVore. Nonlinear approximation. Acta Numerica, 7:51-150, 1998.

[0182] 9. R. A. DeVore and G. G. Lorcutz. Constructive Approximation, Springer-Verlag, New York, 1993.

[0183] 10. D. Donoho. Wedgelets: Nearly-minimax estimation of edges. Annals of Statistics, 27:859-897, 1999.

[0184] 11. D. L. Donoho and X. Huo. Beamlet pyramids: A new form of multiresolution analysis, suited for extracting lines, curves, and objects from very noisy image data. In Proceedings of SPIE 2000, volume 4119, 2000.

[0185] 12. A. C. Gilbert, S. Guha, P. Indyk, Y. Kotidis, S. Muthukrishnan, and M. Strauss. Fast, small-space algorithms for approximate histogram maintenance. In Proceedings of ACM STOC 2002_(—)2002.

[0186] 13. A. C. Gilbert, S. Guha, P. Indyk, S. Muthukrishnan, and M. Strauss. Near-optimal sparse Fourier representations via sampling. In Proc. of ACM STOC, 2002.

[0187] 14. P. Indyk. High-Dimensional Computational Geometry. PhD thesis, Stanford, 2000.

[0188] 15. P. Indyk and R. Motwani. Approximate nearest neighbors: Towards removing the curse of dimensionality. In Proceedings of ACM STOC 1998, pages 604-613, 1998.

[0189] 16. The Journal of Approximation Theory. On the world wide web at math.ohio-state.edu/JAT.

[0190] 17. E. Kushilevitz, R. Ostrovsky, and Y. Rabani. Efficient search for approximate nearest neighbor in high dimensional spaces. In Proceedings of ACM STOC 1998, pages 614-623, 1998.

[0191] 18. S. Mallat and Z. Zhang. Matching pursuits with time frequency dictionaries. IEEE Transactions on Signal Proceeding, 41(12):3397-3415, 1993.

[0192] 19. C. M. Thiele and L. F. Villemoes. A fast algorithm for adapted time frequency tilings. Applied and Computational Harmonic Analysis, 3:91-99, 1996.

[0193] 20. N. Nisan and A. Wigderson. Hardness vs. randomness. J. Comput. System Sci., 49:149-167, 1994.

[0194] 21. Y. C. Pati, R. Rezaiifar, and P. S. Krishnaprasad. Orthogonal Matching Pursuit: Recursive Function Approximation with Applications to Wavelet Decomposition. In Proc. of the 27th Annual Asilomar Conference on Signals, Systems, and Computers, pages 40-44, 1993.

[0195] 22. V. N. Temlyakov. The best m-term approximation and greedy algorithms. Advances in Computational Math., 8:249-265, 1998.

[0196] 23. V. N. Temlyakov. Greedy algorithms and m-term approximation with regard to redundant dictionaries. J. Approximation Theory, to appear.

[0197] 24. Lars Villemoes. Best approximation with Walsh atoms. Constructive Approximation, 13:329-355, 1997.

[0198] 25. Lars Villemoes. Nonlinear approximation with Walsh atoms. In A. Le Mehaute, C. Rabut, and L. L. Schumaker, editors, Surface Fitting and Multiresolution Methods, pages 329-336. Vanderbilt University Press, 1997. 

What is claimed is:
 1. A method comprising a plurality of activities, comprising: automatically: receiving a signal; and via a plurality of iterations, determining, for the received signal, a provable near-optimal B-term representation formed from a linear combination of a plurality of elements selected from a predetermined redundant μ-coherent dictionary, said determined near-optimal B-term representation corresponding to an error that is a minimum of a set of errors corresponding to all possible B-term representations associated with the predetermined redundant μ-coherent dictionary.
 2. The method of claim 1, further comprising: determining B.
 3. The method of claim 1, further comprising: obtaining B.
 4. The method of claim 1, further comprising: determining the redundant μ-coherent dictionary.
 5. The method of claim 1, further comprising: obtaining the redundant μ-coherent dictionary.
 6. The method of claim 1, further comprising: preprocessing the redundant μ-coherent dictionary.
 7. The method of claim 1, further comprising: determining an error of each iterative B-term representation.
 8. The method of claim 1, further comprising: outputting the determined near-optimal B-term representation.
 9. The method of claim 1, wherein the signal is an audio signal.
 10. The method of claim 1, wherein the signal is a speech signal.
 11. The method of claim 1, wherein the signal is a video signal.
 12. The method of claim 1, wherein the near-optimal B-term representation is determined in a time proportional to N raised to a power of approximately 2, where N is a total number of unit vectors comprised by the predetermined redundant μ-coherent dictionary.
 13. The method of claim 1, wherein the near-optimal B-term representation is determined in O(N²) time, where N is a total number of unit vectors comprised by the predetermined redundant μ-coherent dictionary.
 14. The method of claim 1, wherein said determining activity is performed without entirely searching the predetermined redundant μ-coherent dictionary.
 15. The method of claim 1, wherein B is less than 1/(32μ).
 16. The method of claim 1, wherein in O(B/ε) iterations an associated intermediate representation is within an error, denoted 1+E, of the near-optimal B-term representation.
 17. A machine-readable medium comprising instructions for activities comprising: automatically: receiving a signal; and iteratively determining, for the received signal, a provable near-optimal B-term representation formed from a linear combination of a plurality of elements selected from a predetermined redundant μ-coherent dictionary, said determined near-optimal B-term representation corresponding to an error that is a minimum of a set of errors corresponding to all possible B-term representations associated with the predetermined redundant μ-coherent dictionary.
 18. A system comprising: a signal receiver; and a signal processor adapted to automatically iteratively determine, for the received signal, a provable near-optimal B-term representation formed from a linear combination of a plurality of elements selected from a predetermined redundant μ-coherent dictionary, said determined near-optimal B-term representation corresponding to an error that is a minimum of a set of errors corresponding to all possible B-term representations associated with the predetermined redundant μ-coherent dictionary. 